#include <iostream>
#include <Eigen/Dense>

using namespace Eigen;

int main()
{
  MatrixXd m(2,2);
  m(0,0) = 3;
  m(1,0) = 2.5;
  m(0,1) = -1;
  m(1,1) = m(1,0) + m(0,1);
  std::cout << "Here is the matrix m:\n" << m << std::endl;
  VectorXd v(2);
  v(0) = 4;
  v(1) = v(0) - 1;
  std::cout << "Here is the vector v:\n" << v << std::endl;
}

///*
// * Test.cpp
// *
// *  Created on: 23 May 2011
// *      Author: allan
// */
//
//// C++ includes
//#include <cmath>
//#include <iostream>
//#include <iomanip>
//using namespace std;
//
//// FluidReact includes
//#include "AqueousPhase.h"
//#include "System.h"
//#include "Speciation.h"
//#include "Reaction.h"
//#include "BatchReactor.h"
//
//void TestActivityCoefficients_NaCl()
//{
//	AqueousPhase aqueousPhase;
//	
//	aqueousPhase.AddSpecies("H2O(l)");
//	aqueousPhase.AddSpecies("H+");
//	aqueousPhase.AddSpecies("OH-");
//	aqueousPhase.AddSpecies("Na+");
//	aqueousPhase.AddSpecies("Cl-");
//	aqueousPhase.AddSpecies("NaCl(aq)");
//	
//	aqueousPhase.SetElectrolyte("NaCl(aq)", vector<pair<string,double>>({{"Na+", +1.0}, {"Cl-", -1.0}}));
//	
//	aqueousPhase.SetActivityCoefficientModelSetschenow("NaCl(aq)", 0.25);
//	
//	cout << setw(15) << "MolalityNaCl";
//	cout << setw(15) << "IonicStrength";
//	cout << setw(15) << "ActivityH2O";
//	cout << setw(15) << "GammaNaCl+-";
//	cout << setw(15) << "GammaH[+]";
//	cout << setw(15) << "GammaOH[-]";
//	cout << setw(15) << "GammaNa[+]";
//	cout << setw(15) << "GammaCl[-]";
//	cout << setw(15) << "GammaNaCl(aq)" << endl;
//		
//	for(double mNaCl = 0.1; mNaCl < 6.0; mNaCl += 0.1)
//	{
//		aqueousPhase.SetMolesOfSpecies({55.508, 0.0, 0.0, mNaCl, mNaCl, 0.0});
//		
//		vector<double> activityCoeffs = aqueousPhase.ActivityCoefficients(25.0, 100.0);
//		
//		cout << fixed << showpoint;
//		
//		cout << setw(15) << mNaCl;
//		cout << setw(15) << aqueousPhase.GetEffectiveIonicStrength();
//		cout << setw(15) << activityCoeffs[0] * 55.508;
//		cout << setw(15) << sqrt(activityCoeffs[3] * activityCoeffs[4]);
//		cout << setw(15) << activityCoeffs[1];
//		cout << setw(15) << activityCoeffs[2];
//		cout << setw(15) << activityCoeffs[3];
//		cout << setw(15) << activityCoeffs[4];
//		cout << setw(15) << activityCoeffs[5];
//		cout << endl;
//	}
//}
//
//void TestSpeciationWaterNaCl()
//{
//	AqueousPhase aqueousPhase;
//	
//	aqueousPhase.AddSpecies("H2O(l)");
//	aqueousPhase.AddSpecies("H+");
//	aqueousPhase.AddSpecies("OH-");
//	aqueousPhase.AddSpecies("Na+");
//	aqueousPhase.AddSpecies("Cl-");
//	aqueousPhase.AddSpecies("NaCl(aq)");
//	
//	aqueousPhase.SetElectrolyte("NaCl(aq)", vector<pair<string,double>>({{"Na+", +1.0}, {"Cl-", -1.0}}));
//	
//	aqueousPhase.SetActivityCoefficientModelSetschenow("NaCl(aq)", 0.25);
//
//	System system;
//
//	system.SetTemperature(25.0);
//	system.SetPressure(100.0);
//	
//	system.AddPhase(aqueousPhase);
//
//	vector<string> primarySpecies = {"H2O(l)", "H+", "Na+", "Cl-"};
//	vector<string> equilibriumSpecies = {"NaCl(aq)", "OH-"};
//	
//	system.SetPrimarySpecies(primarySpecies);
//	system.SetEquilibriumSpecies(equilibriumSpecies);
//	
//	EquilibriumReaction eqReactionNaClaq;
//	
//	eqReactionNaClaq.SetNameOfReactant("NaCl(aq)");
//	eqReactionNaClaq.SetNameOfProducts(vector<string>({"Na+", "Cl-"}));
//	eqReactionNaClaq.SetStoichiometryOfProducts({1.0, 1.0});
//	eqReactionNaClaq.SetEquilibriumConstant({0.7769}, {25.0});
//	
//	EquilibriumReaction eqReactionOH;
//	
//	eqReactionOH.SetNameOfReactant("OH-");
//	eqReactionOH.SetNameOfProducts(vector<string>({"H2O(l)", "H+"}));
//	eqReactionOH.SetStoichiometryOfProducts({1.0, -1.0});
//	eqReactionOH.SetEquilibriumConstant({13.9951}, {25.0});
//	
//	system.AddEquilibriumReaction(eqReactionNaClaq);
//	system.AddEquilibriumReaction(eqReactionOH);
//	
////	double C_H2O  = 55508;
////	double C_H    = 1.0E-7;
////	double C_OH   = 1.0E-7;
////	double C_Na   = 1.0E-7;
////	double C_Cl   = 1.0E-7;
////	double C_NaCl = 1000;
//	
//	double C_H2O  = 55508;
//	double C_H    = 1.0E-4;
//	double C_OH   = 1.0E-4;
//	double C_Na   = 960;
//	double C_Cl   = 960;
//	double C_NaCl = 40;
//	
//	system.SetMolesOfPhase(0, {C_H2O, C_H, C_OH, C_Na, C_Cl, C_NaCl});
//	
//	Speciation speciation;
//	
//	speciation.SetSystem(system);
//	
//	speciation.ImposeUj(0, 55.508);
//	speciation.ImposeUj(1, 0.0);
//	speciation.ImposeUj(2, 1000);
//	//speciation.ImposeUj(3, U_Cl);
//	
//	speciation.ImposeCe(0, 39);
//	
//	speciation.Solve(1.0E-6, 100);
//	
//	cout << system.GetMolesOfSpecies(0) << endl;
//	cout << system.GetMolesOfSpecies(1) << endl;
//	cout << system.GetMolesOfSpecies(2) << endl;
//	cout << system.GetMolesOfSpecies(3) << endl;
//	cout << system.GetMolesOfSpecies(4) << endl;
//	cout << system.GetMolesOfSpecies(5) << endl;
//	
//	cout << endl;
//	
//	cout << system.GetActivityOfSpecies(0) << endl;
//	cout << system.GetActivityOfSpecies(1) << endl;
//	cout << system.GetActivityOfSpecies(2) << endl;
//	cout << system.GetActivityOfSpecies(3) << endl;
//	cout << system.GetActivityOfSpecies(4) << endl;
//	cout << system.GetActivityOfSpecies(5) << endl;
//}
//
//////void TestSpeciationWaterNaClCaCO3()
//////{
//////	AqueousPhase aqueousPhase;
//////	
//////	aqueousPhase.AddSpecies("H2O(l)");
//////	aqueousPhase.AddSpecies("CaCO3(aq)");
//////	aqueousPhase.AddIon( "H[+]", +1);
//////	aqueousPhase.AddIon("OH[-]", -1);
//////	aqueousPhase.AddIon("Na[+]", +1);
//////	aqueousPhase.AddIon("Cl[-]", -1);
//////	aqueousPhase.AddIon("Ca[2+]", +2);
//////	aqueousPhase.AddIon("HCO3[-]", -1);
//////	aqueousPhase.AddComplex("NaCl(aq)", +1.0, "Na[+]", +1.0, "Cl[-]");
//////	
//////	aqueousPhase.SetActivityCoefficientModelHKF();
//////	aqueousPhase.SetActivityCoefficientModelSetschenow("CaCO3(aq)", 0.1);
//////	aqueousPhase.SetActivityCoefficientModelSetschenow("NaCl(aq)", 0.25);
//////
//////	aqueousPhase.SetTemperature(25.0);
//////	aqueousPhase.SetPressure(100.0);
//////	
//////	MineralPhase calcitePhase;
//////	
//////	calcitePhase.AddSpecies("Calcite");
//////	
//////	EquilibriumReaction eqReactionNaClaq;
//////	
//////	eqReactionNaClaq.SetNameOfReactant("NaCl(aq)");
//////	eqReactionNaClaq.SetNameOfProducts(vector<string>({"Na[+]", "Cl[-]"}));
//////	eqReactionNaClaq.SetStoichiometryOfProducts({1.0, 1.0});
//////	eqReactionNaClaq.SetEquilibriumConstant({0.7769}, {25.0});
//////	
//////	EquilibriumReaction eqReactionOH;
//////	
//////	eqReactionOH.SetNameOfReactant("OH[-]");
//////	eqReactionOH.SetNameOfProducts(vector<string>({"H2O(l)", "H[+]"}));
//////	eqReactionOH.SetStoichiometryOfProducts({1.0, -1.0});
//////	eqReactionOH.SetEquilibriumConstant({13.9951}, {25.0});
//////	
//////	EquilibriumReaction eqReactionCaCO3aq;
//////	
//////	eqReactionCaCO3aq.SetNameOfReactant("CaCO3(aq)");
//////	eqReactionCaCO3aq.SetNameOfProducts(vector<string>({"Ca[2+]", "HCO3[-]", "H[+]"}));
//////	eqReactionCaCO3aq.SetStoichiometryOfProducts({1.0, 1.0, -1.0});
//////	eqReactionCaCO3aq.SetEquilibriumConstant({7.09941}, {25.0});
//////	
//////	System system;
//////
//////	system.AddPhase(aqueousPhase);
//////
//////	vector<string> primarySpecies = {"H2O(l)", "H[+]", "Na[+]", "Cl[-]", "Ca[2+]", "HCO3[-]"};
//////	vector<string> equilibriumSpecies = {"NaCl(aq)", "OH[-]", "CaCO3(aq)"};
//////	
//////	system.SetPrimarySpecies(primarySpecies);
//////	system.SetEquilibriumSpecies(equilibriumSpecies);
//////	
//////	system.AddEquilibriumReaction(eqReactionNaClaq);
//////	system.AddEquilibriumReaction(eqReactionOH);
//////	system.AddEquilibriumReaction(eqReactionCaCO3aq);
//////	
//////	double U_H2O      = 55508.0001;
//////	double U_H        = -0.2033585;
//////	double U_Na       = 1000;
//////	double U_Cl       = 1000;
//////	double U_Ca       = 0.2034585;
//////	double U_HCO3     = 0.2034585;
//////	
//////	double C_H2O      = 55508;
//////	double C_H        = 1.0E-10;
//////	double C_OH       = 1.0E-10;
//////	double C_Na       = 1.0E-10;
//////	double C_Cl       = 1.0E-10;
//////	double C_Ca       = 1.0E-10;
//////	double C_HCO3     = 1.0E-10;
//////	double C_NaClaq   = 1000;
//////	double C_CaCO3aq  = 1.0E-10;
//////	
//////	vector<double> Uj = {U_H2O, U_H, U_Na, U_Cl, U_Ca, U_HCO3};
//////	
//////	vector<double> Cj = {C_H2O, C_H, C_Na, C_Cl, C_Ca, C_HCO3};
//////	
//////	vector<double> Ce = {C_NaClaq, C_OH, C_CaCO3aq};
//////	
//////	vector<double> Ck = {};
//////	
//////	BatchReactor batchReactor;
//////	
//////	batchReactor.SetSystem(system);
//////	batchReactor.SetInitialState(Cj, Ce, Ck);
//////	
//////	batchReactor.SolveBatchNonLinearSystem(Uj, Ck, 1.0E-10, 100);
//////	
//////	cout << scientific;
//////	cout << "C[H2O]      :" << batchReactor.GetCj()[0] << endl;
//////	cout << "C[H]        :" << batchReactor.GetCj()[1] << endl;
//////	cout << "C[Na]       :" << batchReactor.GetCj()[2] << endl;
//////	cout << "C[Cl]       :" << batchReactor.GetCj()[3] << endl;
//////	cout << "C[Ca]       :" << batchReactor.GetCj()[4] << endl;
//////	cout << "C[HCO3]     :" << batchReactor.GetCj()[5] << endl;
//////	cout << "C[NaCl(aq)] :" << batchReactor.GetCe()[0] << endl;
//////	cout << "C[OH]       :" << batchReactor.GetCe()[1] << endl;
//////	cout << "C[CaCO3(aq)]:" << batchReactor.GetCe()[2] << endl;
////////	
////////	cout << endl;
////////	
////////	const AqueousPhase* aqphase = static_cast<const AqueousPhase*>(&batchReactor.GetSystem().GetPhase(0));
////////	
////////	cout << "H2O(l)   :" << aqphase->GetMolarFraction(aqphase->GetIndexOfSpecies("H2O(l)")) << endl;
////////	cout << "H[+]     :" << aqphase->GetMolarFraction(aqphase->GetIndexOfSpecies("H[+]")) << endl;
////////	cout << "OH[-]    :" << aqphase->GetMolarFraction(aqphase->GetIndexOfSpecies("OH[-]")) << endl;
////////	cout << "Na[+]    :" << aqphase->GetMolarFraction(aqphase->GetIndexOfSpecies("Na[+]")) << endl;
////////	cout << "Cl[-]    :" << aqphase->GetMolarFraction(aqphase->GetIndexOfSpecies("Cl[-]")) << endl;
////////	cout << "Ca[2+]   :" << aqphase->GetMolarFraction(aqphase->GetIndexOfSpecies("Ca[2+]")) << endl;
////////	cout << "HCO3[-]  :" << aqphase->GetMolarFraction(aqphase->GetIndexOfSpecies("HCO3[-]")) << endl;
////////	cout << "NaCl(aq) :" << aqphase->GetMolarFraction(aqphase->GetIndexOfSpecies("NaCl(aq)")) << endl;
////////	cout << "CaCO3(aq):" << aqphase->GetMolarFraction(aqphase->GetIndexOfSpecies("CaCO3(aq)")) << endl;
////////	
////////	cout << endl;
////////	
////////	cout << "H2O(l)   :" << aqphase->GetMolarFraction(aqphase->GetIndexOfSpecies("H2O(l)")) * aqphase->ActivityCoefficient(aqphase->GetIndexOfSpecies("H2O(l)")) << endl;
////////	cout << "H[+]     :" << aqphase->GetMolarFraction(aqphase->GetIndexOfSpecies("H[+]")) * aqphase->ActivityCoefficient(aqphase->GetIndexOfSpecies("H[+]")) << endl;
////////	cout << "OH[-]    :" << aqphase->GetMolarFraction(aqphase->GetIndexOfSpecies("OH[-]")) * aqphase->ActivityCoefficient(aqphase->GetIndexOfSpecies("OH[-]")) << endl;
////////	cout << "Na[+]    :" << aqphase->GetMolarFraction(aqphase->GetIndexOfSpecies("Na[+]")) * aqphase->ActivityCoefficient(aqphase->GetIndexOfSpecies("Na[+]")) << endl;
////////	cout << "Cl[-]    :" << aqphase->GetMolarFraction(aqphase->GetIndexOfSpecies("Cl[-]")) * aqphase->ActivityCoefficient(aqphase->GetIndexOfSpecies("Cl[-]")) << endl;
////////	cout << "Ca[2+]   :" << aqphase->GetMolarFraction(aqphase->GetIndexOfSpecies("Ca[2+]")) * aqphase->ActivityCoefficient(aqphase->GetIndexOfSpecies("Ca[2+]")) << endl;
////////	cout << "HCO3[-]  :" << aqphase->GetMolarFraction(aqphase->GetIndexOfSpecies("HCO3[-]")) * aqphase->ActivityCoefficient(aqphase->GetIndexOfSpecies("HCO3[-]")) << endl;
////////	cout << "NaCl(aq) :" << aqphase->GetMolarFraction(aqphase->GetIndexOfSpecies("NaCl(aq)")) * aqphase->ActivityCoefficient(aqphase->GetIndexOfSpecies("NaCl(aq)")) << endl;
////////	cout << "CaCO3(aq):" << aqphase->GetMolarFraction(aqphase->GetIndexOfSpecies("CaCO3(aq)")) * aqphase->ActivityCoefficient(aqphase->GetIndexOfSpecies("CaCO3(aq)")) << endl;
////////	
//////}
//////
//void TestBatchReactorWaterNaClCaCO3()
//{
//	AqueousPhase aqueousPhase;
//	
//	aqueousPhase.AddSpecies("H2O(l)");
//	aqueousPhase.AddSpecies("H+");
//	aqueousPhase.AddSpecies("OH-");
//	aqueousPhase.AddSpecies("Na+");
//	aqueousPhase.AddSpecies("Cl-");
//	aqueousPhase.AddSpecies("Ca++");
//	aqueousPhase.AddSpecies("HCO3-");
//	aqueousPhase.AddSpecies("NaCl(aq)");
//	aqueousPhase.AddSpecies("CaCO3(aq)");
//	
//	aqueousPhase.SetActivityCoefficientModelSetschenow("NaCl(aq)", 0.25);
//	aqueousPhase.SetActivityCoefficientModelSetschenow("CaCO3(aq)", 0.10);
//
//	MineralPhase calcitePhase;
//	
//	calcitePhase.AddSpecies("Calcite");
//	
//	EquilibriumReaction eqReactionOH;
//	
//	eqReactionOH.SetNameOfReactant("OH-");
//	eqReactionOH.SetNameOfProducts(vector<string>({"H2O(l)", "H+"}));
//	eqReactionOH.SetStoichiometryOfProducts({1.0, -1.0});
//	eqReactionOH.SetEquilibriumConstant({13.9951}, {25.0});
//	
//	EquilibriumReaction eqReactionNaClaq;
//	
//	eqReactionNaClaq.SetNameOfReactant("NaCl(aq)");
//	eqReactionNaClaq.SetNameOfProducts(vector<string>({"Na+", "Cl-"}));
//	eqReactionNaClaq.SetStoichiometryOfProducts({1.0, 1.0});
//	eqReactionNaClaq.SetEquilibriumConstant({0.7769}, {25.0});
//	
//	EquilibriumReaction eqReactionCaCO3aq;
//	
//	eqReactionCaCO3aq.SetNameOfReactant("CaCO3(aq)");
//	eqReactionCaCO3aq.SetNameOfProducts(vector<string>({"Ca++", "HCO3-", "H+"}));
//	eqReactionCaCO3aq.SetStoichiometryOfProducts({1.0, 1.0, -1.0});
//	eqReactionCaCO3aq.SetEquilibriumConstant({7.09941}, {25.0});
//	
//	MineralReaction minReactionCalcite;
//	
//	minReactionCalcite.SetNameOfReactant("Calcite");
//	minReactionCalcite.SetNameOfProducts(vector<string>({"Ca++", "HCO3-", "H+"}));
//	minReactionCalcite.SetStoichiometryOfProducts({1.0, 1.0, -1.0});
//	minReactionCalcite.SetEquilibriumConstant({1.84377}, {25.0});
//	minReactionCalcite.SetSpecificSurfaceArea(9.8 * 100.09); // 9.8 cm2/g * 100.09 g/mol
//	
//	MineralReactionMechanism calciteMechanism;
//	
//	calciteMechanism.k  = pow(10.0, -3.48);
//	calciteMechanism.Ea = 35.4;
//	
//	minReactionCalcite.AddMechanisms(calciteMechanism);
//
//	System system;
//
//	system.AddPhase(aqueousPhase);
//	system.AddPhase(calcitePhase);
//
//	vector<string> primarySpecies = {"H2O(l)", "H+", "Na+", "Cl-", "Ca++", "HCO3-"};
//	vector<string> equilibriumSpecies = {"OH-", "NaCl(aq)", "CaCO3(aq)"};
//	vector<string> kineticSpecies = {"Calcite"};
//	
//	system.SetPrimarySpecies(primarySpecies);
//	system.SetEquilibriumSpecies(equilibriumSpecies);
//	system.SetKineticSpecies(kineticSpecies);
//	
//	system.AddEquilibriumReaction(eqReactionOH);
//	system.AddEquilibriumReaction(eqReactionNaClaq);
//	system.AddEquilibriumReaction(eqReactionCaCO3aq);
//	
//	system.AddKineticReaction(minReactionCalcite);
//	
//	double C_H2O      = 55508;
//	double C_H        = 1.0E-7;
//	double C_OH       = 1.0E-4;
//	double C_Na       = 1.0E-4;
//	double C_Cl       = 1.0E-4;
//	double C_Ca       = 1.0E-4;
//	double C_HCO3     = 1.0E-4;
//	double C_NaClaq   = 1000;
//	double C_CaCO3aq  = 1.0E-4;
//	
//	double C_Calcite  = 100;
//	
//	system.SetMolesOfPhase(0, {C_H2O, C_H, C_OH, C_Na, C_Cl, C_Ca, C_HCO3, C_NaClaq, C_CaCO3aq});
//	system.SetMolesOfPhase(1, {C_Calcite});
//	
//	BatchReactor batchReactor;
//	
//	batchReactor.SetSystem(&system);
//	
//	batchReactor.Solve(3, 1.0E-5, 1.0E-3, "BatchReactor[Water-NaCl-Calcite][Simulation 1].dat");
//}
//
//void TestSpeciationWaterCO2()
//{
//	AqueousPhase aqueousPhase;
//	
//	aqueousPhase.AddSpecies("H2O(l)");
//	aqueousPhase.AddSpecies( "H+");
//	aqueousPhase.AddSpecies("OH-");
//	aqueousPhase.AddSpecies("HCO3-");
//	aqueousPhase.AddSpecies("CO2(aq)");
//	
//	GaseousPhase gaseousPhase;
//	
//	gaseousPhase.AddSpecies("H2O(g)");
//	gaseousPhase.AddSpecies("CO2(g)");
//	
//	EquilibriumReaction eqReactionOH;
//	
//	eqReactionOH.SetNameOfReactant("OH-");
//	eqReactionOH.SetNameOfProducts(vector<string>({"H2O(l)", "H+"}));
//	eqReactionOH.SetStoichiometryOfProducts({1.0, -1.0});
//	eqReactionOH.SetEquilibriumConstant({13.9951}, {25.0});
//	
//	EquilibriumReaction eqReactionCO2aq;
//	
//	eqReactionCO2aq.SetNameOfReactant("CO2(aq)");
//	eqReactionCO2aq.SetNameOfProducts(vector<string>({"H2O(l)", "H+", "HCO3-"}));
//	eqReactionCO2aq.SetStoichiometryOfProducts({-1.0, 1.0, 1.0});
//	eqReactionCO2aq.SetEquilibriumConstant({-6.35332}, {25.0});
//	
//	EquilibriumReaction eqReactionH2Og;
//	
//	eqReactionH2Og.SetNameOfReactant("H2O(g)");
//	eqReactionH2Og.SetNameOfProducts(vector<string>({"H2O(l)"}));
//	eqReactionH2Og.SetStoichiometryOfProducts({1.0});
//	eqReactionH2Og.SetEquilibriumConstant({1.49929}, {25.0});
//	
//	EquilibriumReaction eqReactionCO2g;
//	
//	eqReactionCO2g.SetNameOfReactant("CO2(g)");
//	eqReactionCO2g.SetNameOfProducts(vector<string>({"H2O(l)", "H+", "HCO3-"}));
//	eqReactionCO2g.SetStoichiometryOfProducts({-1.0, 1.0, 1.0});
//	eqReactionCO2g.SetEquilibriumConstant({-7.82545}, {25.0});
//	
//	System system;
//
//	system.SetTemperature(25.0);
//	system.SetPressure(100.0);
//	
//	system.AddPhase(aqueousPhase);
//	system.AddPhase(gaseousPhase);
//
//	vector<string> primarySpecies = {"H2O(l)", "H+", "HCO3-"};
//	vector<string> equilibriumSpecies = {"OH-", "CO2(aq)", "H2O(g)", "CO2(g)"};
//	
//	system.SetPrimarySpecies(primarySpecies);
//	system.SetEquilibriumSpecies(equilibriumSpecies);
//	
//	system.AddEquilibriumReaction(eqReactionOH);
//	system.AddEquilibriumReaction(eqReactionCO2aq);
//	system.AddEquilibriumReaction(eqReactionH2Og);
//	system.AddEquilibriumReaction(eqReactionCO2g);
//	
//	double C_H2O      = 55508;
//	double C_H        = 1.0E-7;
//	double C_OH       = 1.0E-4;
//	double C_HCO3     = 1.0E-4;
//	double C_CO2aq    = 1.0E-4;
//	
//	double C_H2Og     = 1.0E-4;
//	double C_CO2g     = 5;
//	
//	system.SetMolesOfPhase(0, {C_H2O, C_H, C_OH, C_HCO3, C_CO2aq});
//	system.SetMolesOfPhase(1, {C_H2Og, C_CO2g});
//	
//	Speciation speciation;
//	
//	speciation.SetSystem(system);
//	
//	speciation.ImposeChargeBalance();
//	speciation.ImposeCj(0, 55.508);
//	speciation.ImposeCe(3, 30);
//	
//	speciation.Solve(1.0E-6, 100);
//	
//	cout << system.GetMolesOfSpecies(0) << endl;
//	cout << system.GetMolesOfSpecies(1) << endl;
//	cout << system.GetMolesOfSpecies(2) << endl;
//	cout << system.GetMolesOfSpecies(3) << endl;
//	cout << system.GetMolesOfSpecies(4) << endl;
//	cout << system.GetMolesOfSpecies(5) << endl;
//	cout << system.GetMolesOfSpecies(6) << endl;
//	
//	cout << endl;
//	
//	cout << system.GetActivityOfSpecies(0) << endl;
//	cout << system.GetActivityOfSpecies(1) << endl;
//	cout << system.GetActivityOfSpecies(2) << endl;
//	cout << system.GetActivityOfSpecies(3) << endl;
//	cout << system.GetActivityOfSpecies(4) << endl;
//	cout << system.GetActivityOfSpecies(5) << endl;
//	cout << system.GetActivityOfSpecies(6) << endl;
//	
//	cout << endl;
//	
//	cout << log10(system.GetEquilibriumReactionQuotient(0)) << endl;
//	cout << log10(system.GetEquilibriumReactionQuotient(1)) << endl;
//	cout << log10(system.GetEquilibriumReactionQuotient(2)) << endl;
//	cout << log10(system.GetEquilibriumReactionQuotient(3)) << endl;
//	
//	cout << endl;
//	
//	cout << "I  = " << system.GetAqueousPhase()->GetEffectiveIonicStrength() << endl;
//	cout << "aw = " << system.GetAqueousPhase()->Activities(25.0, 100.0)[0] << endl;
//	
//	cout << endl;
//	
//	cout << system.GetMolesOfSpecies(1) - system.GetMolesOfSpecies(2) - system.GetMolesOfSpecies(3) << endl;
//}
//
//void TestSpeciationWaterNaClCO2()
//{
//	AqueousPhase aqueousPhase;
//	
//	aqueousPhase.AddSpecies("H2O(l)");
//	aqueousPhase.AddSpecies("H+");
//	aqueousPhase.AddSpecies("OH-");
//	aqueousPhase.AddSpecies("Na+");
//	aqueousPhase.AddSpecies("Cl-");
//	aqueousPhase.AddSpecies("HCO3-");
//	aqueousPhase.AddSpecies("NaCl(aq)");
//	aqueousPhase.AddSpecies("CO2(aq)");
//	
//	GaseousPhase gaseousPhase;
//	
//	gaseousPhase.AddSpecies("H2O(g)");
//	gaseousPhase.AddSpecies("CO2(g)");
//	
//	EquilibriumReaction eqReactionOH;
//	
//	eqReactionOH.SetNameOfReactant("OH-");
//	eqReactionOH.SetNameOfProducts(vector<string>({"H2O(l)", "H+"}));
//	eqReactionOH.SetStoichiometryOfProducts({1.0, -1.0});
//	eqReactionOH.SetEquilibriumConstant({13.9951}, {25.0});
//	
//	EquilibriumReaction eqReactionNaClaq;
//	
//	eqReactionNaClaq.SetNameOfReactant("NaCl(aq)");
//	eqReactionNaClaq.SetNameOfProducts(vector<string>({"Na+", "Cl-"}));
//	eqReactionNaClaq.SetStoichiometryOfProducts({1.0, 1.0});
//	eqReactionNaClaq.SetEquilibriumConstant({0.7769}, {25.0});
//	
//	EquilibriumReaction eqReactionCO2aq;
//	
//	eqReactionCO2aq.SetNameOfReactant("CO2(aq)");
//	eqReactionCO2aq.SetNameOfProducts(vector<string>({"H2O(l)", "H+", "HCO3-"}));
//	eqReactionCO2aq.SetStoichiometryOfProducts({-1.0, 1.0, 1.0});
//	eqReactionCO2aq.SetEquilibriumConstant({-6.35332}, {25.0});
//	
//	EquilibriumReaction eqReactionH2Og;
//	
//	eqReactionH2Og.SetNameOfReactant("H2O(g)");
//	eqReactionH2Og.SetNameOfProducts(vector<string>({"H2O(l)"}));
//	eqReactionH2Og.SetStoichiometryOfProducts({1.0});
//	eqReactionH2Og.SetEquilibriumConstant({1.49929}, {25.0});
//	
//	EquilibriumReaction eqReactionCO2g;
//	
//	eqReactionCO2g.SetNameOfReactant("CO2(g)");
//	eqReactionCO2g.SetNameOfProducts(vector<string>({"H2O(l)", "H+", "HCO3-"}));
//	eqReactionCO2g.SetStoichiometryOfProducts({-1.0, 1.0, 1.0});
//	eqReactionCO2g.SetEquilibriumConstant({-7.82545}, {25.0});
//	
//	System system;
//
//	system.SetTemperature(25.0);
//	system.SetPressure(100.0);
//	
//	system.AddPhase(aqueousPhase);
//	system.AddPhase(gaseousPhase);
//
//	vector<string> primarySpecies = {"H2O(l)", "H+", "Na+", "Cl-", "HCO3-"};
//	vector<string> equilibriumSpecies = {"OH-", "NaCl(aq)", "CO2(aq)", "H2O(g)", "CO2(g)"};
//	
//	system.SetPrimarySpecies(primarySpecies);
//	system.SetEquilibriumSpecies(equilibriumSpecies);
//	
//	system.AddEquilibriumReaction(eqReactionOH);
//	system.AddEquilibriumReaction(eqReactionNaClaq);
//	system.AddEquilibriumReaction(eqReactionCO2aq);
//	system.AddEquilibriumReaction(eqReactionH2Og);
//	system.AddEquilibriumReaction(eqReactionCO2g);
//	
//	double C_H2O      = 55.508;
//	double C_H        = 1.0E-4;
//	double C_OH       = 1.0E-4;
//	double C_Na       = 1.0;
//	double C_Cl       = 1.0;
//	double C_HCO3     = 1.0E-4;
//	double C_NaClaq   = 1.0E-4;
//	double C_CO2aq    = 1.0E-4;
//	
//	double C_H2Og     = 1.0E-4;
//	double C_CO2g     = 30;
//	
//	system.SetMolesOfPhase(0, {C_H2O, C_H, C_OH, C_Na, C_Cl, C_HCO3, C_NaClaq, C_CO2aq});
//	system.SetMolesOfPhase(1, {C_H2Og, C_CO2g});
//	
//	Speciation speciation;
//	
//	speciation.SetSystem(system);
//	
//	speciation.ImposeChargeBalance();
//	speciation.ImposeCj(0, 55.508);
//	speciation.ImposeUj(2, 1);
//	speciation.ImposeUj(3, 1);
//	speciation.ImposeCe(4, 30);
//	
//	speciation.Solve(1.0E-6, 100);
//	
//	cout << system.GetMolesOfSpecies(0) << endl;
//	cout << system.GetMolesOfSpecies(1) << endl;
//	cout << system.GetMolesOfSpecies(2) << endl;
//	cout << system.GetMolesOfSpecies(3) << endl;
//	cout << system.GetMolesOfSpecies(4) << endl;
//	cout << system.GetMolesOfSpecies(5) << endl;
//	cout << system.GetMolesOfSpecies(6) << endl;
//	cout << system.GetMolesOfSpecies(7) << endl;
//	cout << system.GetMolesOfSpecies(8) << endl;
//	cout << system.GetMolesOfSpecies(9) << endl;
//	
//	cout << endl;
//	
//	cout << system.GetActivityOfSpecies(0) << endl;
//	cout << system.GetActivityOfSpecies(1) << endl;
//	cout << system.GetActivityOfSpecies(2) << endl;
//	cout << system.GetActivityOfSpecies(3) << endl;
//	cout << system.GetActivityOfSpecies(4) << endl;
//	cout << system.GetActivityOfSpecies(5) << endl;
//	cout << system.GetActivityOfSpecies(6) << endl;
//	cout << system.GetActivityOfSpecies(7) << endl;
//	cout << system.GetActivityOfSpecies(8) << endl;
//	cout << system.GetActivityOfSpecies(9) << endl;
//	
//	cout << endl;
//	
//	cout << log10(system.GetEquilibriumReactionQuotient(0)) << endl;
//	cout << log10(system.GetEquilibriumReactionQuotient(1)) << endl;
//	cout << log10(system.GetEquilibriumReactionQuotient(2)) << endl;
//	cout << log10(system.GetEquilibriumReactionQuotient(3)) << endl;
//	cout << log10(system.GetEquilibriumReactionQuotient(4)) << endl;
//	
//	cout << endl;
//	
//	cout << "I  = " << system.GetAqueousPhase()->GetEffectiveIonicStrength() << endl;
//	cout << "aw = " << system.GetAqueousPhase()->Activities(25.0, 100.0)[0] << endl;
//	
//	cout << endl;
//}
//
//void TestSpeciationWaterNaClCaCO3CO2()
//{
//	AqueousPhase aqueousPhase;
//	
//	aqueousPhase.AddSpecies("H2O(l)");
//	aqueousPhase.AddSpecies("H+");
//	aqueousPhase.AddSpecies("OH-");
//	aqueousPhase.AddSpecies("Na+");
//	aqueousPhase.AddSpecies("Cl-");
//	aqueousPhase.AddSpecies("Ca++");
//	aqueousPhase.AddSpecies("HCO3-");
//	aqueousPhase.AddSpecies("NaCl(aq)");
//	aqueousPhase.AddSpecies("CO2(aq)");
//	aqueousPhase.AddSpecies("CaCO3(aq)");
//	
//	GaseousPhase gaseousPhase;
//	
//	gaseousPhase.AddSpecies("H2O(g)");
//	gaseousPhase.AddSpecies("CO2(g)");
//	
//	EquilibriumReaction eqReactionOH;
//	
//	eqReactionOH.SetNameOfReactant("OH-");
//	eqReactionOH.SetNameOfProducts(vector<string>({"H2O(l)", "H+"}));
//	eqReactionOH.SetStoichiometryOfProducts({1.0, -1.0});
//	eqReactionOH.SetEquilibriumConstant({13.9951}, {25.0});
//	
//	EquilibriumReaction eqReactionNaClaq;
//	
//	eqReactionNaClaq.SetNameOfReactant("NaCl(aq)");
//	eqReactionNaClaq.SetNameOfProducts(vector<string>({"Na+", "Cl-"}));
//	eqReactionNaClaq.SetStoichiometryOfProducts({1.0, 1.0});
//	eqReactionNaClaq.SetEquilibriumConstant({0.7769}, {25.0});
//	
//	EquilibriumReaction eqReactionCO2aq;
//	
//	eqReactionCO2aq.SetNameOfReactant("CO2(aq)");
//	eqReactionCO2aq.SetNameOfProducts(vector<string>({"H2O(l)", "H+", "HCO3-"}));
//	eqReactionCO2aq.SetStoichiometryOfProducts({-1.0, 1.0, 1.0});
//	eqReactionCO2aq.SetEquilibriumConstant({-6.35332}, {25.0});
//	
//	EquilibriumReaction eqReactionCaCO3aq;
//	
//	eqReactionCaCO3aq.SetNameOfReactant("CaCO3(aq)");
//	eqReactionCaCO3aq.SetNameOfProducts(vector<string>({"Ca++", "HCO3-", "H+"}));
//	eqReactionCaCO3aq.SetStoichiometryOfProducts({1.0, 1.0, -1.0});
//	eqReactionCaCO3aq.SetEquilibriumConstant({7.09941}, {25.0});
//	
//	EquilibriumReaction eqReactionH2Og;
//	
//	eqReactionH2Og.SetNameOfReactant("H2O(g)");
//	eqReactionH2Og.SetNameOfProducts(vector<string>({"H2O(l)"}));
//	eqReactionH2Og.SetStoichiometryOfProducts({1.0});
//	eqReactionH2Og.SetEquilibriumConstant({1.49929}, {25.0});
//	
//	EquilibriumReaction eqReactionCO2g;
//	
//	eqReactionCO2g.SetNameOfReactant("CO2(g)");
//	eqReactionCO2g.SetNameOfProducts(vector<string>({"H2O(l)", "H+", "HCO3-"}));
//	eqReactionCO2g.SetStoichiometryOfProducts({-1.0, 1.0, 1.0});
//	eqReactionCO2g.SetEquilibriumConstant({-7.82545}, {25.0});
//	
//	System system;
//
//	system.SetTemperature(25.0);
//	system.SetPressure(100.0);
//	
//	system.AddPhase(aqueousPhase);
//	system.AddPhase(gaseousPhase);
//
//	vector<string> primarySpecies = {"H2O(l)", "H+", "Na+", "Cl-", "Ca++", "HCO3-"};
//	vector<string> equilibriumSpecies = {"OH-", "NaCl(aq)", "CO2(aq)", "CaCO3(aq)", "H2O(g)", "CO2(g)"};
//	
//	system.SetPrimarySpecies(primarySpecies);
//	system.SetEquilibriumSpecies(equilibriumSpecies);
//	
//	system.AddEquilibriumReaction(eqReactionOH);
//	system.AddEquilibriumReaction(eqReactionNaClaq);
//	system.AddEquilibriumReaction(eqReactionCO2aq);
//	system.AddEquilibriumReaction(eqReactionCaCO3aq);
//	system.AddEquilibriumReaction(eqReactionH2Og);
//	system.AddEquilibriumReaction(eqReactionCO2g);
//	
//	double C_H2O      = 55.508;
//	double C_H        = 1.0E-4;
//	double C_OH       = 1.0E-4;
//	double C_Na       = 1.0E-4;
//	double C_Cl       = 1.0E-4;
//	double C_Ca       = 1.0E-4;
//	double C_HCO3     = 1.0E-4;
//	double C_NaClaq   = 1.0E-4;
//	double C_CO2aq    = 1.0E-4;
//	double C_CaCO3aq  = 1.0E-4;
//	
//	double C_H2Og     = 1.0E-4;
//	double C_CO2g     = 30;
//	
//	system.SetMolesOfPhase(0, {C_H2O, C_H, C_OH, C_Na, C_Cl, C_Ca, C_HCO3, C_NaClaq, C_CO2aq, C_CaCO3aq});
//	system.SetMolesOfPhase(1, {C_H2Og, C_CO2g});
//	
//	Speciation speciation;
//	
//	speciation.SetSystem(system);
//	
//	speciation.ImposeChargeBalance();
//	speciation.ImposeCj(0, 55.508);
//	speciation.ImposeUj(2, 1);
//	speciation.ImposeUj(3, 1);
//	speciation.ImposeUj(4, 1);
//	speciation.ImposeCe(5, 100);
//	
//	speciation.Solve(1.0E-6, 100);
//	
//	cout << setw(10) << system.GetSpecies(0).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(0)  << endl;
//	cout << setw(10) << system.GetSpecies(1).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(1)  << endl;
//	cout << setw(10) << system.GetSpecies(2).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(2)  << endl;
//	cout << setw(10) << system.GetSpecies(3).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(3)  << endl;
//	cout << setw(10) << system.GetSpecies(4).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(4)  << endl;
//	cout << setw(10) << system.GetSpecies(5).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(5)  << endl;
//	cout << setw(10) << system.GetSpecies(6).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(6)  << endl;
//	cout << setw(10) << system.GetSpecies(7).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(7)  << endl;
//	cout << setw(10) << system.GetSpecies(8).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(8)  << endl;
//	cout << setw(10) << system.GetSpecies(9).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(9)  << endl;
//	cout << setw(10) << system.GetSpecies(10).GetName() << ": " << setw(15) << system.GetMolesOfSpecies(10) << endl;
//	cout << setw(10) << system.GetSpecies(11).GetName() << ": " << setw(15) << system.GetMolesOfSpecies(11) << endl;
//	
//	cout << endl;
//	
//	cout << log10(system.GetEquilibriumReactionQuotient(0)) << endl;
//	cout << log10(system.GetEquilibriumReactionQuotient(1)) << endl;
//	cout << log10(system.GetEquilibriumReactionQuotient(2)) << endl;
//	cout << log10(system.GetEquilibriumReactionQuotient(3)) << endl;
//	cout << log10(system.GetEquilibriumReactionQuotient(4)) << endl;
//	cout << log10(system.GetEquilibriumReactionQuotient(5)) << endl;
//}
//
//void TestBatchReactorWaterNaClCaCO3CO2()
//{
//	AqueousPhase aqueousPhase;
//	
//	aqueousPhase.AddSpecies("H2O(l)");
//	aqueousPhase.AddSpecies("H+");
//	aqueousPhase.AddSpecies("OH-");
//	aqueousPhase.AddSpecies("Na+");
//	aqueousPhase.AddSpecies("Cl-");
//	aqueousPhase.AddSpecies("Ca++");
//	aqueousPhase.AddSpecies("HCO3-");
//	aqueousPhase.AddSpecies("NaCl(aq)");
//	aqueousPhase.AddSpecies("CaCO3(aq)");
//	aqueousPhase.AddSpecies("CO2(aq)");
//
//	aqueousPhase.SetActivityCoefficientModelSetschenow("NaCl(aq)", 0.25);
//	aqueousPhase.SetActivityCoefficientModelSetschenow("CaCO3(aq)", 0.10);
//	aqueousPhase.SetActivityCoefficientModelSetschenow("CO2(aq)", 0.20);
//	
//	GaseousPhase gaseousPhase;
//	
//	gaseousPhase.AddSpecies("H2O(g)");
//	gaseousPhase.AddSpecies("CO2(g)");
//	
//	MineralPhase calcitePhase;
//	
//	calcitePhase.AddSpecies("Calcite");
//	
//	EquilibriumReaction eqReactionOH;
//	
//	eqReactionOH.SetNameOfReactant("OH-");
//	eqReactionOH.SetNameOfProducts(vector<string>({"H2O(l)", "H+"}));
//	eqReactionOH.SetStoichiometryOfProducts({1.0, -1.0});
//	eqReactionOH.SetEquilibriumConstant({13.9951}, {25.0});
//	
//	EquilibriumReaction eqReactionNaClaq;
//	
//	eqReactionNaClaq.SetNameOfReactant("NaCl(aq)");
//	eqReactionNaClaq.SetNameOfProducts(vector<string>({"Na+", "Cl-"}));
//	eqReactionNaClaq.SetStoichiometryOfProducts({1.0, 1.0});
//	eqReactionNaClaq.SetEquilibriumConstant({0.7769}, {25.0});
//	
//	EquilibriumReaction eqReactionCaCO3aq;
//	
//	eqReactionCaCO3aq.SetNameOfReactant("CaCO3(aq)");
//	eqReactionCaCO3aq.SetNameOfProducts(vector<string>({"Ca++", "HCO3-", "H+"}));
//	eqReactionCaCO3aq.SetStoichiometryOfProducts({1.0, 1.0, -1.0});
//	eqReactionCaCO3aq.SetEquilibriumConstant({7.09941}, {25.0});
//	
//	EquilibriumReaction eqReactionCO2aq;
//	
//	eqReactionCO2aq.SetNameOfReactant("CO2(aq)");
//	eqReactionCO2aq.SetNameOfProducts(vector<string>({"H2O(l)", "H+", "HCO3-"}));
//	eqReactionCO2aq.SetStoichiometryOfProducts({-1.0, 1.0, 1.0});
//	eqReactionCO2aq.SetEquilibriumConstant({-6.35332}, {25.0});
//	
//	EquilibriumReaction eqReactionH2Og;
//	
//	eqReactionH2Og.SetNameOfReactant("H2O(g)");
//	eqReactionH2Og.SetNameOfProducts(vector<string>({"H2O(l)"}));
//	eqReactionH2Og.SetStoichiometryOfProducts({1.0});
//	eqReactionH2Og.SetEquilibriumConstant({1.49929}, {25.0});
//	
//	EquilibriumReaction eqReactionCO2g;
//	
//	eqReactionCO2g.SetNameOfReactant("CO2(g)");
//	eqReactionCO2g.SetNameOfProducts(vector<string>({"H2O(l)", "H+", "HCO3-"}));
//	eqReactionCO2g.SetStoichiometryOfProducts({-1.0, 1.0, 1.0});
//	eqReactionCO2g.SetEquilibriumConstant({-7.82545}, {25.0});
//	
//	MineralReaction minReactionCalcite;
//	
//	minReactionCalcite.SetNameOfReactant("Calcite");
//	minReactionCalcite.SetNameOfProducts(vector<string>({"Ca++", "HCO3-", "H+"}));
//	minReactionCalcite.SetStoichiometryOfProducts({1.0, 1.0, -1.0});
//	minReactionCalcite.SetEquilibriumConstant({1.84377}, {25.0});
//	minReactionCalcite.SetSpecificSurfaceArea(9.8 * 100.09); // 9.8 cm2/g * 100.09 g/mol
//	
//	MineralReactionMechanism calciteMechanism;
//	
//	calciteMechanism.k  = pow(10.0, -3.48);
//	calciteMechanism.Ea = 35.4;
//	
//	minReactionCalcite.AddMechanisms(calciteMechanism);
//
//	System system;
//
//	system.SetTemperature(25.0);
//	system.SetPressure(100.0);
//	
//	system.AddPhase(aqueousPhase);
//	system.AddPhase(gaseousPhase);
//	system.AddPhase(calcitePhase);
//
//	vector<string> primarySpecies = {"H2O(l)", "H+", "Na+", "Cl-", "Ca++", "HCO3-"};
//	vector<string> equilibriumSpecies = {"OH-", "NaCl(aq)", "CaCO3(aq)", "CO2(aq)", "H2O(g)", "CO2(g)"};
//	vector<string> kineticSpecies = {"Calcite"};
//	
//	system.SetPrimarySpecies(primarySpecies);
//	system.SetEquilibriumSpecies(equilibriumSpecies);
//	system.SetKineticSpecies(kineticSpecies);
//	
//	system.AddEquilibriumReaction(eqReactionOH);
//	system.AddEquilibriumReaction(eqReactionNaClaq);
//	system.AddEquilibriumReaction(eqReactionCaCO3aq);
//	system.AddEquilibriumReaction(eqReactionCO2aq);
//	system.AddEquilibriumReaction(eqReactionH2Og);
//	system.AddEquilibriumReaction(eqReactionCO2g);
//	
//	system.AddKineticReaction(minReactionCalcite);
//	
//	double C_H2O     = 55.508;	
//	double C_H       = 1.0E-4;
//	double C_Na      = 1.0E-4;
//	double C_Cl      = 1.0E-4;
//	double C_Ca      = 1.0E-4;
//	double C_HCO3    = 1.0E-4;
//	
//	double C_OH      = 1.0E-4;
//	double C_NaClaq  = 1.0E-4;
//	double C_CaCO3aq = 1.0E-4;
//	double C_CO2aq   = 1.0E-4;
//	double C_H2Og    = 1.0E-4;
//	double C_CO2g    = 1.0E-4;
//	
//	double C_Calcite  = 2.7;
//	
//	system.SetMolesOfSpecies({C_H2O, C_H, C_OH, C_Na, C_Cl, C_Ca, C_HCO3, C_NaClaq, C_CaCO3aq, C_CO2aq, C_H2Og, C_CO2g, C_Calcite});
//	
//	Speciation speciation;
//	
//	speciation.SetSystem(system);
//	
//	speciation.ImposeChargeBalance();
//	speciation.ImposeCj(0, 55.508);	// H2O
//	speciation.ImposeCj(2, 1);		// Na+
//	speciation.ImposeCj(3, 1);		// Cl-
//	speciation.ImposeCj(4, 1);		// Ca++
//	speciation.ImposeCe(5, 2.0);	// CO2(g) (1 mole = 44 grams)
//	
//	speciation.Solve(1.0E-6, 100);
//	
//	cout << endl;
//	
//	cout << setw(10) << system.GetSpecies(0).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(0)  << endl;
//	cout << setw(10) << system.GetSpecies(1).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(1)  << endl;
//	cout << setw(10) << system.GetSpecies(2).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(2)  << endl;
//	cout << setw(10) << system.GetSpecies(3).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(3)  << endl;
//	cout << setw(10) << system.GetSpecies(4).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(4)  << endl;
//	cout << setw(10) << system.GetSpecies(5).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(5)  << endl;
//	cout << setw(10) << system.GetSpecies(6).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(6)  << endl;
//	cout << setw(10) << system.GetSpecies(7).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(7)  << endl;
//	cout << setw(10) << system.GetSpecies(8).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(8)  << endl;
//	cout << setw(10) << system.GetSpecies(9).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(9)  << endl;
//	cout << setw(10) << system.GetSpecies(10).GetName() << ": " << setw(15) << system.GetMolesOfSpecies(10) << endl;
//	cout << setw(10) << system.GetSpecies(11).GetName() << ": " << setw(15) << system.GetMolesOfSpecies(11) << endl;
//	
//	cout << endl;
//	
//	cout << log10(system.GetEquilibriumReactionQuotient(0)) << endl;
//	cout << log10(system.GetEquilibriumReactionQuotient(1)) << endl;
//	cout << log10(system.GetEquilibriumReactionQuotient(2)) << endl;
//	cout << log10(system.GetEquilibriumReactionQuotient(3)) << endl;
//	cout << log10(system.GetEquilibriumReactionQuotient(4)) << endl;
//	cout << log10(system.GetEquilibriumReactionQuotient(5)) << endl;
//	
//	BatchReactor batchReactor;
//	
//	batchReactor.SetSystem(&system);
//	
//	batchReactor.Solve(5, 1.0E-6, 1.0E-3, "BacthReactor-CO2.dat");
//}
//
//void TestBatchReactorWaterNaClCO2CalciteMagnesite()
//{
//	AqueousPhase aqueousPhase;
//	aqueousPhase.AddSpecies("H2O(l)");
//	aqueousPhase.AddSpecies("H+");
//	aqueousPhase.AddSpecies("OH-");
//	aqueousPhase.AddSpecies("Na+");
//	aqueousPhase.AddSpecies("Cl-");
//	aqueousPhase.AddSpecies("Ca++");
//	aqueousPhase.AddSpecies("Mg++");
//	aqueousPhase.AddSpecies("HCO3-");
//	aqueousPhase.AddSpecies("NaCl(aq)");
//	aqueousPhase.AddSpecies("HCl(aq)");
//	aqueousPhase.AddSpecies("NaOH(aq)");
//	aqueousPhase.AddSpecies("CaCO3(aq)");
//	aqueousPhase.AddSpecies("MgCO3(aq)");
//	aqueousPhase.AddSpecies("NaHCO3(aq)");
//	aqueousPhase.AddSpecies("CO2(aq)");
//	aqueousPhase.AddSpecies("CO3--");
//
//	aqueousPhase.SetActivityCoefficientModelSetschenow("NaCl(aq)", 0.25);
//	
//	GaseousPhase gaseousPhase;
//	gaseousPhase.AddSpecies("H2O(g)");
//	gaseousPhase.AddSpecies("CO2(g)");
//	
//	MineralPhase calcitePhase;
//	calcitePhase.AddSpecies("Calcite");
//	
//	MineralPhase magnesitePhase;
//	magnesitePhase.AddSpecies("Magnesite");
//	
//	EquilibriumReaction eqReactionOH;
//	
//	eqReactionOH.SetNameOfReactant("OH-");
//	eqReactionOH.SetNameOfProducts(vector<string>({"H2O(l)", "H+"}));
//	eqReactionOH.SetStoichiometryOfProducts({1.0, -1.0});
//	eqReactionOH.SetEquilibriumConstant({13.9951}, {25.0});
//	
//	EquilibriumReaction eqReactionNaClaq;
//	
//	eqReactionNaClaq.SetNameOfReactant("NaCl(aq)");
//	eqReactionNaClaq.SetNameOfProducts(vector<string>({"Na+", "Cl-"}));
//	eqReactionNaClaq.SetStoichiometryOfProducts({1.0, 1.0});
//	eqReactionNaClaq.SetEquilibriumConstant({0.7769}, {25.0});
//	
//	EquilibriumReaction eqReactionHClaq;
//	
//	eqReactionHClaq.SetNameOfReactant("HCl(aq)");
//	eqReactionHClaq.SetNameOfProducts(vector<string>({"H+", "Cl-"}));
//	eqReactionHClaq.SetStoichiometryOfProducts({1.0, 1.0});
//	eqReactionHClaq.SetEquilibriumConstant({0.69754}, {25.0});
//	
//	EquilibriumReaction eqReactionNaOHaq;
//	
//	eqReactionNaOHaq.SetNameOfReactant("NaOH(aq)");
//	eqReactionNaOHaq.SetNameOfProducts(vector<string>({"Na+", "H2O(l)", "H+"}));
//	eqReactionNaOHaq.SetStoichiometryOfProducts({1.0, 1.0, -1.0});
//	eqReactionNaOHaq.SetEquilibriumConstant({14.18}, {25.0});
//	
//	EquilibriumReaction eqReactionCaCO3aq;
//	
//	eqReactionCaCO3aq.SetNameOfReactant("CaCO3(aq)");
//	eqReactionCaCO3aq.SetNameOfProducts(vector<string>({"Ca++", "HCO3-", "H+"}));
//	eqReactionCaCO3aq.SetStoichiometryOfProducts({1.0, 1.0, -1.0});
//	eqReactionCaCO3aq.SetEquilibriumConstant({7.09941}, {25.0});
//	
//	EquilibriumReaction eqReactionMgCO3aq;
//	
//	eqReactionMgCO3aq.SetNameOfReactant("MgCO3(aq)");
//	eqReactionMgCO3aq.SetNameOfProducts(vector<string>({"Mg++", "HCO3-", "H+"}));
//	eqReactionMgCO3aq.SetStoichiometryOfProducts({1.0, 1.0, -1.0});
//	eqReactionMgCO3aq.SetEquilibriumConstant({7.34394}, {25.0});
//	
//	EquilibriumReaction eqReactionNaHCO3aq;
//	
//	eqReactionNaHCO3aq.SetNameOfReactant("NaHCO3(aq)");
//	eqReactionNaHCO3aq.SetNameOfProducts(vector<string>({"Na+", "HCO3-"}));
//	eqReactionNaHCO3aq.SetStoichiometryOfProducts({1.0, 1.0});
//	eqReactionNaHCO3aq.SetEquilibriumConstant({0.24474}, {25.0});
//	
//	EquilibriumReaction eqReactionCO2aq;
//	
//	eqReactionCO2aq.SetNameOfReactant("CO2(aq)");
//	eqReactionCO2aq.SetNameOfProducts(vector<string>({"H2O(l)", "H+", "HCO3-"}));
//	eqReactionCO2aq.SetStoichiometryOfProducts({-1.0, 1.0, 1.0});
//	eqReactionCO2aq.SetEquilibriumConstant({-6.35332}, {25.0});
//	
//	EquilibriumReaction eqReactionCO3;
//	
//	eqReactionCO3.SetNameOfReactant("CO3--");
//	eqReactionCO3.SetNameOfProducts(vector<string>({"HCO3-", "H+"}));
//	eqReactionCO3.SetStoichiometryOfProducts({1.0, -1.0});
//	eqReactionCO3.SetEquilibriumConstant({10.32666}, {25.0});
//	
//	EquilibriumReaction eqReactionH2Og;
//	
//	eqReactionH2Og.SetNameOfReactant("H2O(g)");
//	eqReactionH2Og.SetNameOfProducts(vector<string>({"H2O(l)"}));
//	eqReactionH2Og.SetStoichiometryOfProducts({1.0});
//	eqReactionH2Og.SetEquilibriumConstant({1.49929}, {25.0});
//	
//	EquilibriumReaction eqReactionCO2g;
//	
//	eqReactionCO2g.SetNameOfReactant("CO2(g)");
//	eqReactionCO2g.SetNameOfProducts(vector<string>({"H2O(l)", "H+", "HCO3-"}));
//	eqReactionCO2g.SetStoichiometryOfProducts({-1.0, 1.0, 1.0});
//	eqReactionCO2g.SetEquilibriumConstant({-7.82545}, {25.0});
//	
//	MineralReaction minReactionCalcite;
//	
//	minReactionCalcite.SetNameOfReactant("Calcite");
//	minReactionCalcite.SetNameOfProducts(vector<string>({"Ca++", "HCO3-", "H+"}));
//	minReactionCalcite.SetStoichiometryOfProducts({1.0, 1.0, -1.0});
//	minReactionCalcite.SetEquilibriumConstant({1.84377}, {25.0});
//	minReactionCalcite.SetSpecificSurfaceArea(9.8 * 100.09); // 9.8 cm2/g * 100.09 g/mol
//	
//	MineralReactionMechanism calciteMechanism;
//	
//	calciteMechanism.k  = pow(10.0, -5.81);
//	calciteMechanism.Ea = 23.5;
//	
//	minReactionCalcite.AddMechanisms(calciteMechanism);
//	
//	MineralReaction minReactionMagnesite;
//	
//	minReactionMagnesite.SetNameOfReactant("Magnesite");
//	minReactionMagnesite.SetNameOfProducts(vector<string>({"Mg++", "HCO3-", "H+"}));
//	minReactionMagnesite.SetStoichiometryOfProducts({1.0, 1.0, -1.0});
//	minReactionMagnesite.SetEquilibriumConstant({2.03574}, {25.0});
//	minReactionMagnesite.SetSpecificSurfaceArea(9.8 * 84.31); // 9.8 cm2/g * 84.31 g/mol
//	
//	MineralReactionMechanism magnesiteMechanism;
//	
//	magnesiteMechanism.k  = pow(10.0, -9.34);
//	magnesiteMechanism.Ea = 23.5;
//	
//	minReactionMagnesite.AddMechanisms(magnesiteMechanism);
//	
//	System system;
//
//	system.SetTemperature(25.0);
//	system.SetPressure(100.0);
//	
//	system.AddPhase(aqueousPhase);
//	system.AddPhase(gaseousPhase);
//	system.AddPhase(calcitePhase);
//	system.AddPhase(magnesitePhase);
//
//	vector<string> primarySpecies = {"H2O(l)", "H+", "Na+", "Cl-", "Ca++", "Mg++", "HCO3-"};
//	vector<string> equilibriumSpecies = {"OH-", "NaCl(aq)", "HCl(aq)", "NaOH(aq)", "CaCO3(aq)", "MgCO3(aq)", "NaHCO3(aq)", "CO2(aq)", "CO3--", "H2O(g)", "CO2(g)"};
//	vector<string> kineticSpecies = {"Calcite", "Magnesite"};
//	
//	system.SetPrimarySpecies(primarySpecies);
//	system.SetEquilibriumSpecies(equilibriumSpecies);
//	system.SetKineticSpecies(kineticSpecies);
//	
//	system.AddEquilibriumReaction(eqReactionOH);
//	system.AddEquilibriumReaction(eqReactionNaClaq);
//	system.AddEquilibriumReaction(eqReactionHClaq);
//	system.AddEquilibriumReaction(eqReactionNaOHaq);
//	system.AddEquilibriumReaction(eqReactionCaCO3aq);
//	system.AddEquilibriumReaction(eqReactionMgCO3aq);
//	system.AddEquilibriumReaction(eqReactionNaHCO3aq);
//	system.AddEquilibriumReaction(eqReactionCO2aq);
//	system.AddEquilibriumReaction(eqReactionCO3);
//	system.AddEquilibriumReaction(eqReactionH2Og);
//	system.AddEquilibriumReaction(eqReactionCO2g);
//	
//	system.AddKineticReaction(minReactionCalcite);
//	system.AddKineticReaction(minReactionMagnesite);
//	
//	vector<double> C(system.GetNumOfSpecies(), 1.0E-7);
//	
//	C[system.GetIndexOfSpecies("H2O(l)")]    = 55.508;
//	C[system.GetIndexOfSpecies("Na+")]       = 1.0;
//	C[system.GetIndexOfSpecies("Cl-")]       = 1.0;
//	C[system.GetIndexOfSpecies("Ca++")]      = 1.0E-7;
//	C[system.GetIndexOfSpecies("Mg++")]      = 1.0E-7;
//	C[system.GetIndexOfSpecies("HCO3-")]     = 2.0;
//	C[system.GetIndexOfSpecies("CO2(g)")]    = 2.0;
//	C[system.GetIndexOfSpecies("Calcite")]   = 2.0; // Moleweight: 100.09 g/mol, volumic mass: 2709.89 kg/m3
//	C[system.GetIndexOfSpecies("Magnesite")] = 2.0; // Moleweight: 84.31 g/mol, volumic mass: 3009.28 kg/m3
//	
//	system.SetMolesOfSpecies(C);
//	
//	Speciation speciation;
//	
//	speciation.SetSystem(system);
//	
//	speciation.ImposeChargeBalance();
//	speciation.ImposeCj(0, 55.508);	// H2O
//	speciation.ImposeCj(2, 1.0);	// Na+
//	speciation.ImposeCj(3, 1.0);	// Cl-
//	speciation.ImposeCj(4, 1.0E-7);	// Ca++
//	speciation.ImposeCj(5, 1.0E-7);	// Mg++
//	speciation.ImposeCe(10, 2.0);	// CO2(g) (1 mole = 44 grams)
//	
//	speciation.Solve(1.0E-6, 100);
//	
//	cout << endl;
//	
//	for(unsigned iSpecies = 0; iSpecies < system.GetNumOfSpecies(); ++iSpecies)
//		cout << setw(10) << system.GetSpecies(iSpecies).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(iSpecies)  << endl;
//	
//	cout << endl;
//	
//	for(unsigned iEquilibrium = 0; iEquilibrium < system.GetNumOfEquilibriumSpecies(); ++iEquilibrium)
//		cout << "K[" << system.GetSpecies(system.GetIndexOfEquilibriumSpecies(iEquilibrium)).GetName()  << "]: " << setw(15) << log10(system.GetEquilibriumReactionQuotient(iEquilibrium)) << endl;
//	
//	BatchReactor batchReactor;
//	
//	batchReactor.SetSystem(&system);
//	
//	double hour = 3600.0;
//	double day  = 24 * hour;
//	
//	batchReactor.Solve(5 * day, 1.0, 2.0, "BacthReactor.Calcite-Magnesite-CO2-Water-NaCl.dat");
//}
//
//int main()
//{
//	//TestActivityCoefficients_NaCl();
//	
//	//TestSpeciationWaterNaCl();
//	
//	//TestSpeciationWaterNaClCaCO3();
//	
//	//TestBatchReactorWaterNaClCaCO3();
//	
//	//TestBatchReactorWaterNaClCaCO3CO2();
//	
//	TestBatchReactorWaterNaClCO2CalciteMagnesite();
//	
//	//TestSpeciationWaterCO2();
//	
//	//TestSpeciationWaterNaClCO2();
//	
//	//TestSpeciationWaterNaClCaCO3CO2();
//}

